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ABSTRACT 


This thesis develops an iterative algorithm for the design of ARMA models of 
signals in the time domain. ‘The algorithm is based on optimization techniques. 
particularly a gradient technique known as the restricted step method is used. The 
new algorithm is called the iteratzve Prony method, and the results obtained using this 
new method are compared to those obtained using the iterative prefiltering algorithm. 
The thesis shows that the performance of the iterative Prony method is in most of 


the cases comparable or superior to that of the iterative prefiltering algorithm. 


il 


If. 


III. 


TV: 


TABLE OF CONTENTS 


INTRODUCTION .... . . . . SR 
A. THE IDEA OF ARMA MODELING 333 eee 
B. WHY AN ITERATIVE ALGORITHM IN THE SIGNAL DOMAIN 
C. THESIS OUTLINE . .. 23 Se 
ARMA MODELING OF SIGNALS .................. 
A. MODELING METHODS USED IN THIS) tits 
B. OVERVIEW OF MODELING METHODS > eee 

1. Prony’s Method . 2... 

2. Signal Domain Form of Prony s Wetec ree 

3. Iterative Prefiltermg™ Senne 
ITERATIVE ALGORITHM IN THE SIGNAL DOMAIN .... 
A. MULTIDIMENSIONAL OPTIMIZATION BY GRADIENT METH- 


B. THE ITERATIVE PRONY METHOD 22 


1. Vector g oi first derivatives. |. . 7. 592) ee 
2. Matrix G of second derivatives) >. 7) i) ee 
3. Algorithm implementation .... | 92.2332) 
PERFORMANCE OF THE ALGORITHM AND MODELING 
RESULTS ...........2. 32) 335 ee 
A. TEST DATA USED IN THIS THESIS 3 
B. PRESENTATION OF RESULTS) ee 


1. Simulated test data -...... . . [ee 
2 Acoustic test data... «<0 6 eco 5 
CONCLUSIONS ..............,, 3 


Pees CUSSION OF RESULTS 


B. RECOMMENDATIONS FOR FUTURE WORK 


LIST OF REFERENCES ..... 


Ce eC eT Ce ee ee | 


. os « e) ses “e s&s «ss cs os eos o8 68 8 oe oe e¢ oe so oe «oe «6 


e ee ee @e@ o#  @e@ ee ®  ®#  @®  @®  e# ee oe oe @e  e# @e¢ ee «se ee 


LIST OF TABLES 


DESCRIPTION OF SIMULATED TES 28 
DESCRIPTION OF ACOUSTIC TEST DATA 
POLE-ZERO LOCATION OF THE SYSTEMS USED TO MODEL 
THE SIMULATED NOISY TEST DATA 


V1 


Zt 
Zee. 
2.3 
2.4 


Sat 
3.2 


3.3 


4.1 


4.2 


LIST OF FIGURES 


Block diagram for the direct method for signal modeling. ....... 
Block diagram for the indirect modeling problem. ........... 
Block diagram for the iterative prefiltering method. .......... 
Signal wren01 and its 4 poles/4 zeros models. (a) Using Prony’s 
method. (b) Using Signal Domain form of Prony’s method. (c) Using 


Prerahive a MCMUCHING =e ne ee 


Displacement of the poles and zeros of an iterative Prony’s model 

Displacement of the poles and zeros of a 4-3 model. The modeled 
signal in this case actually has two poles on the real axis ....... 
Displacement of the poles and zeros of a 4-3 model. The initial model 
shows two poles in the real axis, the final model in this case has no 


POlCoOMMIMGNCHOMGGEeE Fu Ges UR dee doe ee 


Signal t01_n and its 2 poles-1 zero models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
tO1_n and the iterative prefiltering model. (c) Signal t01_n and the 
What ChimEOIVaTOdeMmr sj 5 40.5 oe ee ek 
signal t01_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
t01_n and the iterative prefiltering model. (c) Signal t01_n and the 


itchoni cumnonvmodeimayeee ny OEP SGSS Se el ee. 


Vu 


Cyt 


CO 


Lt 


13 


20 


4.3 


4.4 


4.6 


4.7 


Signal t02_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
(02_n and the iterative prefiltering model. (c) Signal t02_n and the 
iterative Prony model. ... ... 2 ee 
Signal ¢02_n and its 6 poles-5 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
t02_n and the iterative prefiltering model. (c) Signal t02_n and the 
iterative Prony model... 4 i952 
Signal t03_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
t03_n and the iterative prefiltering model. (c) Signal t03_n and the 
iterative Prony model. .. . . . . ge 
Signal t03_n and its 6 poles-5 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
t03_n and the iterative prefiltering model. (c) Signal t03_n and the 
iterative Prony model. .. . . 2. . 1 3 epee neerennen 
Signal t04_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
tO4_n and the iterative prefiltering model. (c) Signal t04_n and the 
iterative Prony model. ..... . . . (ee 
Signal t04_n and its 6 poles-5 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
tO04_n and the iterative prefiltering model. (c) Signal t04_n and the 


iterative Prony model. ...... . . . sey 


Vill 


4.9 


4.10 


4.11 


Ael2 


4.13 


4.14 


Signal ¢05_n and its 2 poles-1 zero models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
t05.n and the iterative prefiltering model. (c) Signal t05_n and the 
wach Pivorty ane), 9 es, 
Signal t05_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
t05_n and the iterative prefiltering model. (c) Signal t05_n and the 
iterahiveE ion y Inede a cee. 2... 
Signal vowel_a being modeled by iterative prefiltering using a (10,9) 
order system. (a) Normalized squared-norm of the error between the 
model and the actual signal. (b) Signal vowel_a and the iterative 
prefiltering model selected from the 20‘ iteration. (c) Signal vowel_a 
and the iterative prefiltering model selected from the 17‘ iteration. 

Behavior of the poles of the system used to model the signal vowel_a 
during one oscillation of the iterative prefiltering algorithm. (a) 15*” 
iteration. (b) 16%” iteration. (c) 17" iteration. (d) 18" iteration. (e) 
emieranon 20 aiteraiion em a ye ge 
Behavior of the zeros of the system used to model the signal vowel_a 
during one oscillation of the iterative prefiltering algorithm. (a) 15% 
iteration. (b) 16 iteration. (c) 17 iteration. (d) 18" iteration. (e) 


eo aeiccra Tony 20 MacratiOn ee eKeee. . . 


Signal wren0/ and its 7 poles-6 zeros models. (a) Normalized squared- 


norm of the error between the models and the actual signal. (b) Signal 
wren0l and the iterative prefiltering model. (c) Signal wren@1 and the 


ey CulanommyeiMOdelme ses fg. 2 a 


1X 


43 


4.19 


4.20 


Signal wren@1 and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual sig- 
nal. (b) Signal wren@/ and the iterative prefiltering model. (c) Signal 
wren01 and the iterative Prony model. (2). e ee 
Signal vowel_a and its 10 poles-9 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. 
(b) Signal vowel_a and the iterative prefiltering model. (c) Signal 
vowel_a and the iterative Prony model. (=) 3 eee 
Signal vowel_a and its 14 poles-13 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. 
(b) Signal vowel_a and the iterative prefiltering model. (c) Signal 
vowel_a and the iterative Prony model... yyy see 
Signal b20_2133a and its 8 poles-7 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. 
(b) Signal bi0_2133a and the iterative prefiltering model. (c) Signal 
b1o-2133a and the iterative Prony model. 7.5). see ee 
Signal b20_2133a and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. 
(b) Signal bi0_2133a and the iterative prefiltering model. (c) Signal 
bi0_2733a and the iterative Prony model, oa.s =e 
Signal bio_2385a and its 8 poles-7 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. 
(b) Signal b70_2385a and the iterative prefiltering model. (c) Signal 


bi0o_2385a and the iterative Prony model. 2 9 s)aueuene een 


30 


ol 


a2] 


4,23 


Signal bt0_2385a and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. 
(b) Signal bio0_2385a and the iterative prefiltering model. (c) Signal 
pions es omamanune Itenaibivese bony IMOdels.5.. . 646 5 se we el. 
Signal 60_80a and its 8 poles-7 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal 
bio_80a and the iterative prefiltering model. (c) Signal bio_80a and the 
Rei PATO MCCS a a 
Signal bz0_80a and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual sig- 
nal. (b) Signal bio_80a and the iterative prefiltering model. (c) Signal 


Messvamimauheiterative Prony modelo... . 1 et ee. 


X] 


ACKNOWLEDGMENT 


To my children, Carlos Andres, Diego Felipe , and Maria Carolina whose con- 
stant efforts to prevent me from finishing this thesis fortunately failed. To my wife, 
Monica who deserves most of the credit for this work. Only her love, dedication. 
inspiration, and support during these last two years made all this possible. I am also 
indebted to my thesis advisor, Dr. Charles W. Therrien, for introducing me to the 
fields of signal processing and for pointing me in the right direction in the completion 


of this work. 


XU 


I. INTRODUCTION 


A. THE IDEA OF ARMA MODELING 


The goal of linear modeling is to accurately represent an observed data sequence 
as the output of a linear filter. The idea of representing a complicated process with a 
comparatively simpler model has many different applications. Curve fitting in math- 
ematical modeling, analysis of electronic devices using equivalent circuits, and system 
transfer functions in automatic control are just a few examples outside the area of dig- 
ital signal processing. Parametric modeling has also a large number of applications in 
signal processing. Currently there is considerable interest in the parametric modeling 
approach to spectral estimation. In speech processing the applications include digital 
transmission, storage, and synthesis of the speech signal. Our particular interest is 
the modeling of sonar signals, such as biologics and other underwater acoustic data. 
This work forms part of an overall research program in sonar signal modeling. The 
research will help to understand the relative benefits of signal domain algorithms 
versus algorithms based on coefficients of the transfer function. It is hoped that the 
method developed in this thesis will become an important tool in the overall effort 
for sonar signal modeling. 

In linear modeling the filter used to generate the data sequence is usually rep- 
resented by a linear difference equation with constant coefficients. The z-transform 
of this type of system is a rational polynomial function. Three type of models are 
derived from this kind of systems; they are known as autoregresive (AR), moving 
average (MA), and autoregresive moving average (ARMA). 

Much work has been done on AR models, which correspond to all-pole systems. 


The reason for that is that the parameters for the model can be obtained by solving 


linear equations, and a great body of theory has been developed that applies to 
this problem [{Ref. 1, 2]. Relatively much less work has been done with MA and 
ARMA models. However. since MA models have limited applications and are almost 
as difficult to obtain as ARMA models, most interest centers on the latter. The 
filter in this type of models has both poles and Zeros, and the design fundamentally 
involves nonlinear equations. A properly designed ARMA model can provide better 
performance than an AR model, with a smaller number of parameters. ARMA 


modeling is the topic of this thesis. 


B. WHY AN ITERATIVE ALGORITHM IN THE SIG- 
NAL DOMAIN 


There are many different approaches to the problem of ARMA modeling. The 
majority of them are based on statistical techniques [Ref. 3]. Some of these methods 
regard the data as a realization of a random process while others focus on the data as 
given (Ref. 1]. Data oriented methods try to minimize some criterion that estimates 
how well the model fits the data, in most cases the least squares error between the 
signal and the model. Stochastic approaches may attempt to estimate the model 
parameters directly from the data by solving nonlinear equations or by spectral fac- 
torization. The maximum likelihood procedure, for example [Ref. 1, 4], is essentially 
nonlinear. A number of indirect methods have been developed that modify the norm 
of the error by separating the AR and MA parts of the problem so that at least some 
of the equations to estimate the parameters become linear. This approach is found 
in procedures such as the Prony’s method, Shank’s method, and the least squares 
modified Yule-Walker method [Ref. 1, 2, 5, 6]. A different type of approach replaces 
the nonlinear problem with iteration while trying to solve for the AR and MA pa- 
rameters simultaneously. The iterative prefiltering method of Steiglitz and McBride 


(Ref. 7, 8] is of this type. 


Our method is also of the latter type. However, the advantage is that it works 
directly with the poles of the rational model. which we know affect the performance 
of the system. The poles are displaced in specific directions so that the new model 
minimizes the error between the model output and the original signal. Iterative 
prefiltering on the other hand works with the coefficients of the transfer function, 
so it is difficult or impossible to predict its effects on the poles and zeros of the 
system. The new algorithm is much more dependable with respect to convergence 
than the iterative prefiltering algorithm because it moves poles and zeros specifically 


to minimize the error between the model and the original signal. 


C. THESIS OUTLINE 


The remainder of this thesis is organized as follows. Chapter II presents the 
modeling methods that are used in this thesis and gives a brief explanation of all 
of them. Chapter III introduces the reader to the theory of multidimensional op- 
timization by gradient methods and develops the iterative Prony method, which is 
the main contribution of this thesis. Chapter IV presents the results of testing the 
algorithm on simulated and real acoustic data and compares these results with those 
obtained using iterative prefiltering. Chapter V gives conclusions and suggestions for 


future research. 


II. ARMA MODELING OF SIGNALS 


A. MODELING METHODS USED IN THIS THESIS 


This thesis deals with deterministic approaches to ARMA modeling. Two types 
of modeling methods are considered. First we have non-iterative methods like Prony’s 
method and its alternate signal domain form |Ref. 1]; second we consider iterative 
methods like iterative prefiltering and the new zterative Prony method developed in 
this thesis. 

The goal of Prony’s method is to represent a given sequence z{n] as the impulse 
response of a linear time invariant (LTI) system. In the transform (z) domain this 


representation has the form 


(2.1) 





where X(z) is the z-transform of z[n] and B(z)/A(z) represents the transfer function 
of the system. This approximation is explained in more detail in the next section. 
What has become known as Prony’s method in the current signal processing literature 
differs from Prony’s original work in some respects. Our basic form of Prony’s method 
solves for the coefficients of the transfer function in (2.1). 

The signal domain form of Prony’s method is closer to Prony’s original work 
[Ref. 1, p.560] and seeks to represent the data in terms of a set of damped exponen- 
tials as 


r[n] © crt + cers +--+ + cprp 


where the r; are the roots of the denominator polynomial A(z) and the cx are the 
complex coefficients required for the expansion. Both forms involve linear equations 


and least squares techniques. 









<= ona b,z* 
A(z) 1+ Sp 44274 








Figure 2.1: Block diagram for the direct method for signal modeling. 


The iterative prefiltering method used is due to Steiglitz and McBride [Ref. 7, 8] 
and attempts to match the data with a model of finite order using an efficient iterative 
approach. It differs from Prony’s method in that it solves for both numerator and 
denominator polynomial coefficients stmultaneously at each iteration. 

Whenever a signal is modeled using a fixed-order rational polynomial model, an 
approximation has to be made and some kind of measure has to be used to determine 
the “goodness” of the model. In this thesis the least squares error norm (mbozl, 
norm) is used to measure the approximation error. This norm measures the energy 
of the error and is the norm most widely used primarily due to its mathematical 


tractability [Ref. 2]. 


B. OVERVIEW OF MODELING METHODS 


1. Prony’s Method 


The derivations of all the modeling methods presented in this section follow 
those in [Ref. 1, pp. 550-564] and [Ref. 2]. As stated above, Prony’s method aims 
at representing a signal as the impulse response of an LTI system. Figure 2.1 shows 
an implementation of this approximation which is known as the direct method. The 
system function X(z) in the transform domain is a rational polynomial function 
B(z)/A(z) with Q zeros and P poles. The error signal e[n] is computed as the 


difference between the response of the system z[n] and the given signal z[n], i.e. 


20] Se pl Sa eeltole (2.2) 


The LTI system is chosen to minimize the sum of squared errors 


Sp = >_ |e[n]|’. (2338 


This problem leads to nonlinear equations whose solution (if a unique solution actu- 
ally exists) turns out to be a very difficult task [Ref. 1. 2]. To avoid these difficulties. 
a number of indirect methods have been developed for modeling. Prony’s method is 
one of these procedures and can be derived as follows. The LTI system satisfies the 
difference equation 

t[n] + ayz[n — 1] +--+ + apa[n — P] = bo d{n] + b,6[n — 1] +--- + bgd[n — Q]. 


(2.4) 
If the requirement that 
tin] = cll; 7 Oe ee 


is applied to (2.4), where N, is the length of the data, and the difference equation is 





evaluated for n = 0,1,...,N, — 1, the result is the matrix equation 

x0] 0 0 --- 0 bo 

cl r|0] 0 +++ 0 by 

ee) eal r|0| en 1 by 

‘at 

: : : ee ag : 

; ; | _|: (95 
o1Q) -e(Q-1) 2(Q-2] -- efQ-P| 3 ig | 
Q+1) 21Q]  2lQ-1) ~~ clqg-P+i |]; 0 

: meet ap 

z[N,—1] z{[N,-—2] 2z[N,-3] --- z[Ns-—P—1] | 0 


This can be written as 


aad (2.6) 


where 


z[0] 0 0 0 
zl] 2/0] 0 0 
PA eal] [0] 0 
sci es ae | (2.7) 


LiGheni@ethuai@~ 2) <-> 20 


and 


z[(Q+1) 2|Q] z[(Q-—1) --- 2|Q-P+]] 
X,= | | | —_ (2.8) 
z[N,—-1] z[N,-—2] 2z[N,—-3] --- 2[Ns—P-1] 
and a and b are the vectors of coefficients appearing in (2.5). The lower partition 


Xa ="0 (2.9) 


represents an overdetermined set of linear equations that need not have an exact 
solution. This set of equations can be solved by least squares, where a is chosen to 


minimize the least squares norm of the equation error S4 = |le,||* in 
Xa = C4. (2.10) 


This leads to a set of linear equations called the normal equations, which can be 


written compactly as [Ref. 1, pp.536-537] 
(X7X4)a= | SA | (ain) 


and can be solved for a and S4. Once the vector a is known, it can be substituted 
in the upper partition of (2.6) 


to solve for the vector b. Although it is referred here simply as Prony’s method. 


this procedure is also known as the modern Prony method or the extended Prony 


method. 


ealn] = x[n] * a[n] — d[n] 





bo, bi, «250 bo. OLOnO ne 


Figure 2.2: Block diagram for the indirect modeling problem. 


Although Prony’s method is simple to implement. it is important to keep 
in mind that it is an indirect method. In particular this procedure (see Fig. 2.2) 


minimizes the squared magnitude of the error 
e4(n] = z[n] * a[n| — d[n] @a) 


where the sequences a[n] and b/n] are defined as 


def | @,; 0S See F 
a[n] = 0; otherwise (2.14) 
and 
def } Gn; OS n= Oop — ae 
ae = 0; otherwise (2-19) 


where P and Q are the number of poles and zeros respectively. Equation 2.13 repre- 
sents a different least squares error from that in (2.3) where the quantity to minimize 
was the squared magnitude of the error e|n] (see Fig. 2.1). The practical significance 
of this difference is that frequently there is a loss of accuracy in estimating the poles 


and zeros by Prony’s method [Ref. 3]. 
2. Signal Domain Form of Prony’s Method 


An alternative formulation of the method described above can be obtained 
if the problem (2.1) is stated in the signal domain by representing z|n] in terms of a 


set of complex exponentials: 


t[n] S qr} + carp +--+ + cpr (21a 


where as stated before, the r, are the roots of the polynomial A(z), assumed to be 
distinct, and the complex coefficients c, provide for the linear combination of the P 
roots. 

This approximation can be initially formulated as in the section above and 
(2.9) can be solved in the least squares sense for the coefficients of A(z) (ie. the 
vector a). The roots r; of A(z) can then be found and (2.16) can be evaluated for 


n=0,1,...,N; —1 to produce the set of equations 


C) 
ry r2 bee a C lel 
re r2 oe pt, a = § 7/2) (2.17) 
pNel pNel .,, pNewl | oP | z[Ns — 1] 


This set of equations can then be solved in a least squares sense to obtain the vector 
of coefficients c. 

In the case of multiple roots at the same location, a slight variation of the 
same procedure can be used. Suppose, for example, that r; is a double root. In this 


case, the approximation is 
rin] & ary tenrt +--+ cprp (23) 


and the matrix equation to solve for the coefficients becomes 


1) tee tr x(0] 
Cy 
inl re Ts ..s Pp ; [1] 
ie or he oe Te, ee Tr) (2710) 
ee (N, — rie cae tes rps SP r{Ns — 1] 


This situation is rare, however, because computational errors and errors inherent to 
the modeling method itself contribute to produce roots that may be very close to 


each other, but not exactly at the same location. 


3. Iterative Prefiltering 


The iterative prefiltering method attempts to solve the “direct” problem 
mentioned in subsection 1 of this chapter and to refine the initial pole-zero estimate 
by solving a succession of linear problems. Equation 2.2 (error for the direct problem) 
can be written in the z-domain as 


B(z) _ A(Z)Alz) =e 
A(z) A(z) | 





E(z) = X(z) —-AX@ == (2.20) 


The notion of iteration can be introduced that allows computation of a new set of 
poles and zeros based on the last known set of poles. Iterative prefiltering replaces 
the error for the direct problem at the (2 + 1) iteration with the iterative error 


function 
X(z)A@*Y(z) — BENQ) 
Al) (z) 


If h(n] is the inverse z-transform of 1/A™(z), then the error can be written in the 


ECM (z) = (2.21) 


signal domain as 
et Din} = a[n] * h(n] * at) [n] — b°*) [n] « AM In). (2.228 


The coefficients at!) {n] and b@+){n] are selected to minimize the corresponding sum 


of squared errors 
Ns-l . 9 
SUrU= SP ec aria (2:28) 
mie" 
at each iteration; this situation is shown in Fig. 2.3. No general proof of convergence 
has been given for this algorithm; however, it is easy to see that if the iteration does 
converge, it must produce the same answer as the direct method. Specifically, at 
convergence A‘) = A(t) and (2.21) becomes the same as (2.20). 
If we use an indirect modeling procedure like Prony’s method to compute 


the initial vector a and we define x) as 
cin} & [nj * h[n), (2.24) 
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Figure 2.3: Block diagram for the iterative prefiltering method. 


then the sequences h{[n] and x“ [n] for n = 0,1,..., Ns — 1 can be computed from 


the recursive difference equations 
h(n] = 6[n] — sofa [n — k] (2.25) 
and 
z'[n) = z[n] — > oe 2) [n — k]. (2,26) 


Thus the error (2.22) can be written as 


Q . : 
elt In 2 a) 2O[n — k] — So oS MAO In — ji]. ON 
j7=0 
In order to find the coefficients a! +1) and aes the error is written in matrix form 


ior 2 = 0) 1,...,NV,=— las 


2 a ali+}) a 
xO HO | | cs | (41) (2.28) 
where ; 
eM a) Pa 2)/P — 1) x)(0] 
x() e9Y[P+1] 2 [P] cI (2.29) 
eINs—1) x [Ng — 2] 2) {Ns —1— P] 
h®|P] AOlP -— 1 hO(P —Q] 
wo | 4 P+) AMP) + AYP -Q+]] (2.30) 
hOINs — 1) hO[LNs — 2] hYINs —1—Q] 


and 


1 pert) 
(t+1) (7+1) 
q(t) = ‘. ) b(t!) — 1 . (238 
i+1) i+] 
als ae ) 


Equation 2.28 is analogous to (2.10). Thus the least squares problem defined by 


minimizing the norm of the error in (2.28) can be reduced to 


S+1) 
: xT ; ; a(t) 0 
([ xX Ae | xO FAL 1) | pen | = |" (2.32) 
0 
where 
Sot) = eae (2.33) 


This linear equation can then be solved for the vectors of filter parameters. 


At this point, it is instructive to compare the performance of the three 
modeling methods outlined in this chapter by applying each to model a small segment 
of a transient sound corresponding to a wrench being struck. This wrench sound was 
recorded and sampled in the laboratory at a sampling rate of approximately 10,240 
Hz. This signal, denoted by wren01, was also used and modeled in [Ref. 9]. Figure 
2.4 shows a 100 point segment of the signal wren01 with one of the three models 
(Prony’s, signal domain form of Prony’s, and Iterative Prefiltering ) overlaid. In all 
three cases the signal was modeled with four poles and four zeros. As can be seen, 
the difference between iterative prefiltering and the first two models is significant. 
The non-iterative methods match only the initial points of the sequence, and produce 
poor approximations in modeling the remaining part of the signal. This is due, in 
large part, to poles not sufficiently close to the unit circle. Iterative prefiltering, 
on the other hand, produces a model which is close to the real sequence along the 


complete segment. 






a aaa “A 
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Figure 2.4: Signal wren0Q1 and its 4 poles/4 zeros models. (a) Using Prony’s method. 
(b) Using Signal Domain form of Prony’s method. (c) Using Iterative Prefiltering. 
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Ill. ITERATIVE ALGORITHM IN THE 
SIGNAL DOMAIN 


A. MULTIDIMENSIONAL OPTIMIZATION BY GRADI- 
ENT METHODS 


A multidimensional function f(z), 22,-::,2,) that is continuous and differen- 
tiable can be minimized using one of several very powerful hillclimbing techniques 
known as gradient methods [Ref. 10, p. 84]. Some of those methods are derived on 
the basis of a quadratic model that can be obtained from a truncated Taylor series 
expansion of f(x). Let x‘*) denote the value of x at the k'" iteration. Then for any 


point x = x*) + 6; when 6 is small, the function can be approximated by 
| 1 
f(x +5) x g(6) =f +g" 5 + =67GM5 (3.1) 


where g and G represent the vector of first derivatives and the matrix of second 
derivatives of the function f(x) respectively and they should be available at every 
point. In Newton’s method the iterate x(*+)) is taken to be x(*) + 6“), where the 
correction 5‘) minimizes q‘*)(5). This method is only well defined when the matrix 
of second derivatives G is positive definite, in which case the k“" iteration of Newton’s 
method is given by the following procedure [Ref. 11, pp 44-46]: 

1. solve G5 = —g() for 6 = 5") 

2. set xfFt1) = xf) 4 5 (4), 


(3.2) 


The fact that G() may not be positive definite when x“) is far from the solution, 
and that even when G) is positive definite convergence may not occur, makes this 


method undesirable as a general formulation of a minimization algorithm. However, a 
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number of variations to the basic method have been proposed that are more suitable 
for a general class of problems. One of these methods is Newton's method with line 
search [Ref. 11, pp. 47-49] in which the Newton algorithm is used to generate a 


direction of search 


o(k) — Bet) ot*) (3.3) 


which can later be used in a line search algorithm to actually calculate the correction 
§. In the cases when G) is not positive definite, the linear search can be made 
along +s‘*) choosing the correct sign to ensure a descent direction. However, some 
difficulties that arise here (like very high numerical costs and failure of convergence 
for some special cases) make this an undesirable approach for our algorithm. 

As stated before, Newton’s method is defined only when the matrix G *) is 
positive definite, and this matrix is positive definite only when the error 6 is “small”: 
or better stated, the method is defined only in some neighborhood 2) of x) in 
which q‘*)(§) agrees with f(x) + 6) in some sense. In such cases, it is correct to 
choose x(*+!) = x(4) 4 §) with the correction 6“) minimizing g‘*)(6) for all x) +6 
in QQ"), This method is referred to as the restricted step method because the step is 
restricted by the region of validity of the Taylor series. [Ref. 11] 


The region of definition for the k** iteration can be expressed as 


Gi) Spe 2 een (3.4) 
where || - || denotes the norm of the vector. In this case, the optimization problem 
can be stated as: 

minimize, q')(6) subject to |||] < h®. (3.5) 


As mentioned before, the least squares norm is the one most commonly used in this 
type of problems, so it is the one used in this thesis and is denoted as || - ||). The 


problem that now becomes apparent is how to select the error margin h'*) of the 


1d 


neighborhood (3.4). This margin should be as large as possible because the iteration 
step is directly related to it. Various methods have been proposed to control the 
parameter h'*); one of these methods attempts to insure that the Newton’s search 
direction problem (3.3) is always defined [Ref. 11. pp. 100-103]. It does so by adding 


a multiple of the unit matrix I to G“) and computing the new problem 
(G® + rT) 6 = —g (3.6) 


where the net effect is that increases in y cause ||6||2 to decrease, and vice versa. 


If we define 





ay SFO _ flO — fort + 6) 
Ag) ft) — g(64)) 


(3.7) 
then the ratio r*) represents a measure of the accuracy to which q'*)(6*) approxi- 
mates f(x) + 6)) on the k** step, and as the accuracy increases r‘*) gets closer to 
unity. Using (3.7), Marquardt [Ref. 12] suggests an algorithm that tries to adaptively 
maintain h(*) as large as possible while controlling the ratio r“). The k'® iteration 
of such an algorithm is stated as: 

1. given x) and v*), calculate g“) and G™; 

2. factor G&) + yI; if not positive definite. reset v“) = 4v'*) and repeat; 

3. solve (3.6) to find 5“): 


4. evaluate f(x) + 6)) and hence r; 


5. if rf) < 0.25 set vt!) = 4p) 
else if r®) > 0.75 set vt) = y\*) /9 
else set v(4+)) = (4). 


6. if r(*) < Q set xl4*t)) = x(*) else set x(4t) = x(4) 4 6), 


Here the parameters 0.25, 0.75, 4 and 2 are arbitrary, and v™) > 0 is also chosen 
arbitrarily [Ref. 11, pp. 102-103]. Proofs of global and second order convergence for 
this algorithm are given in [Ref. 11, pp. 96-98] for the cases when the first and sec- 


ond derivatives of the function f(x) exist, and the vector x“) belongs to a bounded 
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n-dimensional space for all &. Although this method does have some disadvantages. 
it represents a good basis for the formulation of a general minimization algorithm. 
Some variations of this method were considered, but it was found that for the spe- 
cific application of ARMA modeling in the time domain, these variations produced 


an extremely high overhead in calculations. 
B. THE ITERATIVE PRONY METHOD 


Let us now return to the problem of representing a sequence z[n] as a linear 


combination of complex exponentials. Equation 2.17 can be written as 
Re=x-+e (3.8) 


where e is called the equation error, x represents the data which may or may not 
be complex, c is the vector of complex coefhcients, and R is the matrix of complex 


roots, which can be written more specifically as 


lk ] l 
TR, + Jia Ge ITI ei) ea ITIp 
AP re) aos (TRp ale hates) (3.9) 


i (rR, ar jrn)? (Gee 


N.-1 N,-1 


Mgr tn) * (ta, - rr) (rrp + grip) | 


where rp, and ry represent the real and imaginary components of the 2‘ root, re- 
spectively, and N, is the number of data samples. 


By defining 
OQ |lelly = «77 = (Rc —x)*” (Re — x), (3.10) 


it is clear that the problem is to find the vector r= ry+ jr; of P complex roots that 


minimizes the function Q(r). 


The first and second derivatives of Q with respect to r are represented by the 
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vector g and matrix G.which are defined by 

















39 
Orr, 
ag 
OrR, 
90 
—_ OTR 
S=\| ag (3.11) 
Orr 
Ig 
Or 
a0 
OrIp 
and 
_8QO #4 —89Q |... _#O 7 oo 2. ee 
Orr, Orr, Orr, OTR, Orr, OTRp Or R, Ort, “ae iE OrR,OTIp 
Orr, Orr, Orr, OTR, OTR, OTR p OrR, OT), ECT Orr, rrp 
ae 079 — a0 070 a0 — 070 
Orp Orp Orr OrR OrR OTR OrR Ory Orp OTT, OrR Ory 
= P i P 2 J gd bee 1 2 P Ie 
G=| #9" #07" | #9" #9) #0" ., _#e! 
Orr, Orr, Orr, OrR, or, OTRp One Orr, Ory, Ory, arr, Ort, OTT» 
_#9' _#9° | _#g° _#@° _#9° | _ sg” 
Ory,Orr, Or, OTR, Shape Ta Or, Orr, a Ory, Orrp 
_#oQ _@9Q |... _ #9 £—_§00 oo a 
L Orrporr, OTT, OTR, OrrpOTRp Orr, Ort, OrTp Ors, Orr pOTtp 
(3.12) 


In order to provide for a more compact notation, define the gradient operator with 


respect to the complex vector r as the vector of partial derivatives 











Vee} JF | =| ee (3.13) 
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consequently, (3.11) can be expressed as 




























































































Ves) 
| iis / 
a = Vi, — : (3.14) 
Vr, 2 
Equation 3.12 can also be written as 
ra) JO dO JQ ake IQ ex?) 
Orr, Orr, OrR, OTRp Ory, Or I5 OTT p 
a dQ dQ ex® AQ JQ e} 
Orr, Orr, Orr, OTRp shale Or I, Ont» 
3 a9 82 «., 82 8s QQ a9 | 
OrRp Orr, OTR, OTRp ory, on Orr p 
a= (ley 
a exe) dQ dQ Ie) dQ, exe) 
Orr, Orr, 9IFTR, Orrp Ory, Ori OFT 
3 dQ dO 3Q exe) 3Q a0 | 
Ort, Orr, or Rp ee ary, OTT, Orrp 
e) dQ dQ dQ Jo 00 39 | 
Orr p Orr, OIFTR, Orrep Orr, OF Or tp 








o_ | 8 al 








Orr Orr, Orr, 
G= 5 
a. /{ 22 2g 
OF I, Orr, ony, 
Gee a 
bel (3.16) 


and it becomes clear that the matrix of second derivatives can be expressed compactly 


as 


eo eta, 2), (Vr,2) |. (3.17) 


The following subsections derive explicit expressions for the two quantities g and G. 
1. Vector g of first derivatives 


From (3.8) it is easy to see that 


«TT «T 
us me d cm eee (3.18) 


Or; 7 Or; ca Or; Or; 


ly) 


where r; represents the real or imaginary parts of the 7" root. Using (3.18) and the 


chain rule in (3.10) leads to 


dQ ., OR Ones 
Orn, =e€E Orr, c+ Cc Orr, Ex (3, 19) 








and explicit evaluation of the partial derivatives of the matrix R results in 





0 
I 
AQ i 2(rR, + ITY, ) 
Orr, — 3 (rR, Sa “ + 


(N, =1) (ie, + gee mee 


a | 0 1 2(rr,—Jrr) 3(rr,—Jrn)? --+ (Ns—V (rr, — arn)? | ; 


(3.20) 


where c¢; is the 2" component of the vector c. If the quantity ¢; is defined as 











0 
l 
2(rr, + Jr) 
def ' ' 
=a 3(rp, +irn) os (3.21) 
NU ) j Ns—2 
(N; —1) (rr, +3r1) 
then 
OR 
— sen ‘ 
Orr Cc S12) (3 22) 
ONRE 
ew = ¢77, (3.23) 
Orr, 
and (3.20) can be written compactly as 
0 
= eT ¢. ae (3.24) 
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At this point it can be recognized that (3.24) represents the addition of two scalar 
quantities where the first one is the complex conjugate of the second, so (3.24) can 


be further simplified to 





rn = 2 Re |¢i7e| (3.25) 


where Re[-] denotes the real part of the vector. Finally, using (3.25) for 2 = 1---P 
in the upper partition of (3.14) produces 
on 
gy 
VreY = Z Re E€ . (3.26) 
ip 
A similar procedure can be used to obtain the gradient of Q with respect 


to the imaginary part of the vector r. Once more, from (3.18) and the chain rule 


applied to (3.10), it follows that 


OR , 
OR* 
«T BI sh *T 
Cc ten ye; (3.28) 


These results can be used to generate an expression similar to that in (3.24) for the 


vector of first derivatives of Q with respect to the imaginary components of r 


JQ __ -,x«x] 
Ben Be. Glee (3.29) 


which in turn defines the vector of partial derivatives with respect to the imaginary 


components as 

T 
G1 
«T 
oo 


Vr,2 =2Im eS. (3.30) 


«xT 
Gp 


Equations 3.26 and 3.30 can now be combined as shown in (3.14) to obtain the final 


expression for the vector of first derivatives 


Re E 


(3.31) 


2. Matrix G of second derivatives 


An expression for the matrix G of second derivatives can be obtained as 


follows. Substituting (3.24) and (3.29) into (3.16) yields 


wee ( [eer Fei] Flees E77 |) 





sie ( [eTeterte] jfeTe~ei%] ) 


Then using the chain rule and both expressions in (3.18) leads to 


C= 








eT oT 
T O8, 4 oxTORT, 4 ¢+T_aR 38. “| aT 26, ‘TORT. _ .sT OR 7 eee 
E OT Re mee OTR, - a g Or Ry a OTR, : pe Orr, aS Or R,. - G: Orr, 7 Orr, : 











eT oT 
mae 3g, Thorne. xT OR 3g. eee 36, «TOR! . — ,*) ORae 3g, 
c Ort, ae Orr, cit 8 orn, ~ 7 Ort, “ ah Orr, a Orr, oie Orr, © Orr, 1 


(3.33) 


to 
Se) 


From the definition of € in (3.21) it 1s seen that 


0 
0 
2 

Grae c¢;, if i=k 
ne 





=s, = : Ono 
Orr, Sik . | ( ) 


| iN ae 1L)(N, — 2) (rR, ey a | 








0, otherwise. 
In the same way it can be shown that 

aA 

Orr — ote fares) ) 
k 

OE , | 

5 = JSik (3.36) 
k 

oer" - 

a a jie [onan 
k 


These last four expressions together with (3.22), (3.23),(3.27), and (3.28) can now be 
substituted in (3.33) to obtain 


le" eu + ei7e; + ei7e,+eiFe] 5 [e*Tsu + ei"e; -— €77e, — siFe| 


C= 
j le? s, — e476; + €77€, — sie] — [ese — CFP; — 7, + iP 
ol eee 
k=1...P (3.38) 


Finally notice again that the elements of the matrix G are formed by addi- 


tions and subtractions of complex scalars with their respective complex conjugates; 
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thus the final expression for G is given by 
2 Re é77e,] +2Re|sife] 2 Im [e77e,] — 2 Im [sje] 
—2 Im ler7e,| =p Ihae site| Eee ere, | 2 Re [sife| 
(=) eS eee (3.39) 
where it should be remembered that s;, = 0 for all: 4 &. 


3. Algorithm implementation 


An iterative method for ARMA modeling in the time domain was imple- 
mented using the results of the last two subsections in conjunction with the algo- 
rithm presented in section A of this chapter. We call this method the zteratzve Prony 
method. 

The algorithm uses the signal domain form of Prony’s method to calculate 
an initial model for the given sequence. From there it uses the calculated model 
to compute the error e, the vector of first derivatives g, and the matrix of second 
derivatives G and iterates until specific conditions are met. Figure 3.1 is an example 
of how the algorithm changes the position of the poles and zeros of the initial model 
in order to minimize the error. This figure represents the poles and zeros of a transfer 
function of order (4,3)—4 poles 3 zeros—that was overmodeled using a (6,5) order 
model. It is clear from the figure that the tendency in this case is to have a pole-zero 
cancellation (see second and third quadrants) of two poles and two zeros as expected. 

Some features were added to the basic algorithm in order to deal with special 
modeling cases. Specifically, if the initial model has some roots on the real axis, then 
because of the way the algorithm iterates, those roots never move away from the real 
axis. A modification was therefore introduced to deliberately displace those roots 


from the real axis and proceed with the iterations. If the tendency of the roots is 
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Figure 3.1: Displacement of the poles and zeros of an iterative Prony’s model 


“to go back” to the real axis, then they are returned to their initial position, and the 
iterations continue. Otherwise the roots may continue to spread apart and move as 
a complex pair. Figure 3.2 is an example of the displacement of the poles and zeros 
of an order (4,3) model. In this case the modeled signal actually has two poles on 
the real axis, and the initial model correctly placed two of the poles in the real axis. 
Those poles are displaced from the real axis by the algorithm, but then after some 
iterations it is clear that the poles are tending to return to the real axis. At this point 
the poles are forced back to the real axis by setting their imaginary parts to zero 
and the iterations continue until convergence is obtained. The opposite situation is 
shown in Figure 3.3. In this case the initial model also has two poles located on the 
real axis, but contrary to the case presented above, the roots, after being displaced 
from the real axis, continue to move away from the axis until they reach their final 


position in the first and fourth quadrants closer to the unit circle. 


20 


08 * = poles before iterations 
0.6 x = poles during iterations 


0.4 oO = zeros during iterations 


Im{] 


0.2 





-0.4 


Figure 3.2: Displacement of the poles and zeros of a 4-3 model. The modeled signal 
in this case actually has two poles on the real axis 


* = poles before iterations 
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Figure 3.3: Displacement of the poles and zeros of a 4-3 model. The initial model 
shows two poles in the real axis, the final model in this case has no poles in the axis 


IV. PERFORMANCE OF THE ALGORITHM 
AND MODELING RESULTS 


A. TEST DATA USED IN THIS THESIS 


Two types of signals were used in this thesis to test the performance of the 
algorithm. The first type, which will be called szmulated test data, consists of five 
sequences (t01 to t05) each one hundred points long that were produced as the 
impulse response of a known rational system. Noise to produce an SNR in the range 
of 10 to 15 dB was added to the original sequences, and the resulting sequences were 
designated as t01_n to t05.n. The original signals are described in Table 4.1 by their 
transfer functions and the location of their poles and zeros. 

The second group of test signals consists of recorded acoustic data. ‘Two of these 
signals were recorded and sampled in the laboratory. One of them is the sequence 
wren01 already mentioned in Chapter II; the other one was obtained from human 
speech, in particular, the signal vowel_a corresponds to 100 samples of the spanish 
vowel a. The remaining three signals from the group of acoustic data were recorded 
at sea by a submarine platform; they correspond to sounds produced by marine life 


and ice cracking. The description of the acoustic signals is presented in Table 4.2. 
B. PRESENTATION OF RESULTS 


All simulated test signals were modeled twice with the iterative Prony method. 
In the first test the exact number of poles and zeros of the original model was used; 
in the second test all signals were modeled using two more poles and zeros than the 
original model. This last test is considered closer to a real life situation where the 


exact order of the signal to be modeled is unknown. 


taf) 


TABLE 4.1: DESCRIPTION OF SIMULATED TEST DATA 
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TABLE 4.2: DESCRIPTION OF ACOUSTIC TEs TDA 


SIGNAL | DESCRIPTION 
Transient sound corresponding to a wrench being struck 
Spanish vowel a 


bio_2133a | Sperm whale 


bio_2385a | Porpoise whistle 


bio_80a Ice cracking 





To mathematically produce a meaningful measure of the performance of a mod- 
eling algorithm can be quite difficult since various norms can be deceiving when 
comparing errors of signals with large differences in magnitude. Two different ap- 
proaches to measure the performance of the algorithms were therefore used in this 
thesis. The first is quantitative and involves computing the squared-norm of the error 
between the model and the actual signal and dividing it by the total energy (norm) 
of the signal. The second approach is to overlay in a plot the model and the original 
signal in order to provide a visual comparison of the results. This is less quantitative 


but frequently more revealing of errors in the modeling process. 
1. Simulated test data 


The first data sets modeled were the simulated test data sets. Figure 4.1(a) 
is a comparison of the normalized errors that result when the sequence t01_n is 
modeled with 2 poles and 1 zero using both iterative prefiltering and iterative Prony 
methods. Figure 4.1(b) and (c) show 100 points of the sequence t01_n and the 
two order (2,1) models. At this point there is no noticeable difference between the 
iterative prefiltering and the iterative Prony models. Figure 4.2(a) again shows a 
comparison of the normalized errors between an iterative prefiltering model and an 
iterative Prony model of t01_n for the case when the signal t01_n was overmodeled 
using models of order (4,3). Although the difference between the two modeled signals 
in this case is not large, notice that the error for the iterative prefiltering method 
initially increases before decreasing while the error for the iterative Prony method 
decreases monotonically. This is the first example of a pattern that repeats in all but 
one of the simulated test signals that were modeled. Figures 4.3 through 4.10 give 
similar comparisons for the remaining simulated signals. The pattern, which can be 
seen in Figures 4.1 to 4.10, is that when the signals are modeled with a number of poles 


and zeros different from that of the actual order of the signal (always overmodeling 
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Figure 4.1: Signal t0/_n and its 2 poles-1 zero models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t0/_n and 
the iterative prefiltering model. (c) Signal {0/_n and the iterative Prony model. 
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Figure 4.2: Signal ¢0/_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t0/_n and 
the iterative prefiltering model. (c) Signal t0i_n and the iterative Prony model. 
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Figure 4.3: Signal {02_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal #02_n and 
the iterative prefiltering model. (c) Signal t02_n and the iterative Prony model. 
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Figure 4.4: Signal ¢02_n and its 6 poles-5 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t02_n and 
the iterative prefiltering model. (c) Signal t02_n and the iterative Prony model. 
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Figure 4.5: Signal t03_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t03_n and 
the iterative prefiltering model. (c) Signal ¢03_n and the iterative Prony model. 
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Figure 4.6: Signal t03_n and its 6 poles-5 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t03_n and 
the iterative prefiltering model. (c) Signal t03_n and the iterative Prony model. 


35 


wee solid = Iterative Prony method 
dashed = Iterative prefiltering method 


jee) Wee 8 Pl 





S 12k ! 
sa cE 5 
= 5 
FO ee SS 
- H 
; 
1033 | | Iteration number | 
0 ] y) 3 4 5 6 
(a) 
2 ) 
15 solid = signal t04_n 
v , dashed = iterative prefiltering model | 
1; 
= \ | 
= 
= 0. | ; 
& 





| n 
10 20 30 40 50 60 70 80 90 100 
(b) 
2 
is solid = signal t04_n 
2 dashed = iterative Prony model 
= 1 
& acl 
2 °° | \ 
0 : ye NN 
a Vas ‘s | 
-0.5 
0 10 20 30 40 50 60 70 80 90 100 
(c) 


Figure 4.7: Signal t04_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal (04_n and 
the iterative prefiltering model. (c) Signal #04_n and the iterative Prony model. 
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Figure 4.8: Signal ¢04_n and its 6 poles-5 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t0/_n and 
the iterative prefiltering model. (c) Signal t04_n and the iterative Prony model. 
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Figure 4.9: Signal ¢05_n and its 2 poles-1 zero models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t05_n and 
the iterative prefiltering model. (c) Signal ¢05_n and the iterative Prony model. 
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Figure 4.10: Signal t05_n and its 4 poles-3 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal t05_n and 
the iterative prefiltering model. (c) Signal t05_n and the iterative Prony model. 
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in this case), the behavior of the iterative prefiltering method tends to degrade (1.e., 
it takes longer for the algorithm to reach convergence) while the behavior of the 
iterative Prony method remains the same or even improves as in the case of t05n 
shown in Figures 4.9(a) and 4.10(a). Parts (b) and (c) of Figures 4.1 through 4.10 
show the original signals and their respective models overlaid. It can be seen that 
the models arrived at by both methods follow the original signals very closely in all 
cases. Table 4.3 lists the location of the poles and zeros of the systems used to model 
all the simulated noisy signals. These systems were obtained using both the iterative 
prefiltering and the iterative Prony methods. The poles and zeros shown in Table 4.3 
can be compared to the poles and zeros of the original simulated signals presented in 
Table 4.1. It is clear that the location of the poles and zeros of the modeled signals 
should not be exactly the same as those of the original signals because some noise (in 
the order of 10 to 15 dB SNR) was intentionally added before the modeling process. 
However, Tables 4.1 and 4.3 show a close relation between the location of the poles 
and zeros of the original sequences and the position of the poles and zeros of the 


modeled signals. 
2. Acoustic test data 


As mentioned in the last section the acoustic test data represents sounds 
recorded both underwater and in a laboratory environment. In some cases shorter 
segments were selected for modeling due to the complexity of these signals. Once 
again the iterative prefiltering algorithm was used, and its results were compared to 
the results obtained using the iterative Prony method. 

Before presentation of results, it is important to explain how the model 
produced by the iterative prefiltering method was selected. For all the cases, the 
same number of iterations was used both for the iterative Prony and for the iterative 


prefiltering methods. In order to obtain the best possible results from the iterative 
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TABLE 4.3: POLE-ZERO LOCATION OF THE SYSTEMS USED TO MODEL 
THE SIMULATED NOISY TEST DATA 


SIGNAL JPR IE a ITERATIVE 
IP 1GCOAN NG PREFILTERING 












NAME & 
(ORDER) | POLES = |ZEROS = |POLES | EZEROST 


Pei ee Sea, | ea ORT oT 
(2,1) 
t01_n 0.9510 Z + 0.6290 | 0.7636 0.9506 Z + 0.6289 | 0.7630 
Prezzo rrr 
t02_n 0.9545 Z + 0.6238 | 13.2102 0.9250 Z + 0.6313 | 6.8483 
| ia) [oisiss coors | 21565, a7 [asian ce 0.0088 | resto: 0.1650 _| 
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prefiltering algorithm (especially in those cases when convergence was not obtained), 
the model selected was that corresponding to the iteration with the lowest error 
between the model and the original sequence. Figure 4.11 is an example of one of 
the cases when convergence was not reached using the iterative prefiltering method. 
It can be seen that if the system selected is the one obtained at the last iteration 
(20'") of the algorithm (which in this case corresponds to a peak of the error), then 
the resulting model does not follows the original signal at all. However, if we select 
the system from the iteration for which the error is the smallest (17 iteration in this 
case), then we obtain a model that is closer to the original signal. Some insight can 
also be obtained if we look at the behavior of the poles and zeros of the system while 
the error of the iterative prefiltering algorithm is oscillating. Figures 4.12 and 4.13 
show the position of the poles and zeros during the oscillation that exists between 
iterations 15 to 20 of the case presented in Figure 4.11. For the first part of these 
iterations when the error is low, the poles remain at approximately the same position. 
At iteration 20, (Fig. 4.12(f)) the poles suddenly spread apart, some to even outside 
the unit circle, causing the large increase in the error. A similar effect occurs with 
the zeros, although there is some significant movement of the zeros in the earlier 
iterations. This type of behavior is avoided in the iterative Prony method where a 
controlled and systematic displacement of the poles and zeros is produced to finally 
reduce the error between the modeled and original signals. 

The same approach used in the last subsection was used to model the acous- 
tic data. All signals were first modeled with a low order system and subsequently 
remodeled using a higher order system. The results are shown in Figures 4.14 to 4.23. 
In most of the cases the models obtained using the iterative Prony method closely 
follow the original signals. It can also be observed that the iterative Prony method 
does not have as large a dependency in the order of the model selected as does the 


iterative prefiltering method. While it is of course necessary to select a model with a 
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Vigure 4.11: Signal vowel_a being modeled by iterative prefiltering using a (10,9) 
order system. (a) Normalized squared-norm of the error between the model and the 
actual signal. (b) Signal vowel_a and the iterative prefiltering model selected from 
the 20%" iteration. (c) Signal vowel_a and the iterative prefiltering model selected 
from the 17 iteration. 
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Figure 4.12: Behavior of the poles of the system used to model the signal vowela 
during one oscillation of the iterative prefiltering algorithm. (a) 15“ iteration. (b) 
16%" iteration. (c) 17** iteration. (d) 18** iteration. (e) 19% iteration. (f) 20” 
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Figure 4.13: Behavior of the zeros of the system used to model the signal vowel_a 
during one oscillation of the iterative prefiltering algorithm. (a) 15% iteration. (b) 
16%" iteration. (c) 17 iteration. (d) 18%" iteration. (e) 19%" iteration. (f) 20% 
iteration. 
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high enough order to properly model the original signal, once such a model has been 
selected, increments or small variations to its order do not produce degradations on 
the performance of the algorithm. On the contrary, it was found that in some cases 
the performance of the iterative prefiltering algorithm can be significantly reduced 
when the order of the model is increased slightly. It can also be seen that the itera- 
tive Prony algorithm tends to provide a closer match to the data than the iterative 
prefiltering method, and in most of the cases the rate of convergence of the iterative 
Prony method was higher than that of iterative prefiltering. Another important point 
that can be extracted from the results presented in Figures 4.14 to 4.23 is that while 
convergence with neither algorithm is guaranteed, we obtained convergence with the 
iterative Prony method in all cases for these acoustic signals. The same was not 
true for iterative prefiltering. In most of the cases the error for the iterative Prony 
method begins to decrease starting at the first iteration, and although the change is 
not monotonic in all cases, the error after a few iterations is consistently lower than 


the initial error. 
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Figure 4.14: Signal wren0/ and its 7 poles-6 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal wren0/ and 
the iterative prefiltering model. (c) Signal wren0/ and the iterative Prony model. 
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Figure 4.15: Signal wren@J and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
wren01 and the iterative prefiltering model. (c) Signal wren@/ and the iterative 
Prony model. 
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Figure 4.16: Signal vowel_a and its 10 poles-9 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
vowel_a and the iterative prefiltering model. (c) Signal vowel_a and the iterative 
Prony model. 
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Figure 4.17: Signal vowel_a and its 14 poles-13 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
vowel_a and the iterative prefiltering model. (c) Signal vowel_a and the iterative 
Prony model. 
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Figure 4.18: Signal 6:0_2133a and its 8 poles-7 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
bio_2133a and the iterative prefiltering model. (c) Signal 620_2/33a and the iterative 
Prony model. 
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Figure 4.19: Signal bto_2133a and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
bi0_2133a and the iterative prefiltering model. (c) Signal bio_2133a and the iterative 
Prony model. 
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Figure 4.20: Signal bi0_2385a and its 8 poles-7 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
bio_2385a and the iterative prefiltering model. (c) Signal b¢o_2385a and the iterative 
Prony model. 
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Figure 4.21: Signal bi0_2385a and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
b10.2385a and the iterative prefiltering model. (c) Signal bi0_2385a and the iterative 
Prony model. 
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Vigure 4.22: Signal bio_80a and its 8 poles-7 zeros models. (a) Normalized squared- 
norm of the error between the models and the actual signal. (b) Signal bo_80a and 
the iterative prefiltering model. (c) Signal bzo_80a and the iterative Prony model. 
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Figure 4.23: Signal bzo_80a and its 12 poles-11 zeros models. (a) Normalized 
squared-norm of the error between the models and the actual signal. (b) Signal 
bio_S0a and the iterative prefiltering model. (c) Signal bz0_80a and the iterative 
Prony model. 
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V. CONCLUSIONS 


A. DISCUSSION OF RESULTS 


In this thesis, a new method for modeling signals in the time domain is developed 
and applied to model both recorded acoustic data and simulated signals produced 
as the impulse response of a known system. We call this method the iterative Prony 
method. In most of the simulated test data sets the models provided by the iterative 
Prony method are sufficiently close to the original signals; in most cases. it is diff- 
cult to distinguish between the signal and the model. When modeling the acoustic 
data distortion becomes apparent in some of the models, which may be due to the 
complexity of the structure of the signals. However, this distortion is no worse than 
for any of the best algorithms that have been used to model this data previously. 

The new algorithm was compared very specifically to the iterative prefiltering 
algorithm (Ref. 7, 8] which has been used in modeling a variety of acoustic data 
[Ref. 3, 9]. The rate of convergence of the iterative Prony method was in most of the 
cases comparable or superior to that of the iterative prefiltering algorithm. Thus, 
while iterative prefiltering sometimes has convergence problems, the new algorithm 
is much more dependable in that respect. The price to pay for this improvement is 
in the number of computations. While the number of floating point operations per 


iteration in the iterative prefiltering method is approximately 
64(P ++Q —1)° + 8N,(P + Q)? + 10(P + Q)N, + 12N,, 
iterative Prony requires about 


672P° + (24N, + 102) P? + (60.N, + 46)P + 198N, 
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floating point operations at each iteration. For example for a complex data set of 
length 100 (NV, = 100) and P = Q = 8 we have approximately 452,400 floating 
point operations per iteration using iterative prefiltering versus 572,360 using iter- 
ative Prony. If we increase the order of the model to P = Q = 16, then we have 
approximately 2,787,824 operations per iteration in the iterative prefiltering algo- 


rithm versus 3,509,560 in the iterative Prony method. 
B. RECOMMENDATIONS FOR FUTURE WORK 


The iterative prefiltering algorithm has been the main tool in the modeling 
efforts for sonar data modeling [Ref. 9, 13]. The new iterative Prony algorithm is 
now at a stage where it can be substituted for the iterative prefiltering algorithm and 
tested in operational use. To do so needs some further programming to make the 
segmentation of the data automatic and to make the entire modeling procedure more 
of a “turn crank” operation. These should be some of the very next steps. In addition, 
the practical implications of the increased computation needs to be addressed, and 
if possible new methods need to be developed to help reduce computations. 

In a larger sense the work reported in this thesis can be used as a base for 
possible applications of the iterative Prony method in the problems of filter design, 
speech processing, and spectral estimation. The expressions for the vector of first 
derivatives g and the matrix of second derivatives G of the error derived in Chapter 
III can be used along with different minimization methods to provide for other new 


modeling methods that may adapt better to specific modeling problems. 
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